Transcriptome analysis of human cholangiocytes exposed to carcinogenic 1,2-dichloropropane in the presence of macrophages in vitro

1,2-Dichloropropane (1,2-DCP), a synthetic organic solvent, has been implicated in causality of cholangiocarcinoma (bile duct cancer). 1,2-DCP-induced occupational cholangiocarcinoma show a different carcinogenic process compared to common cholangiocarcinoma, but its mechanism remains elusive. We reported previously that exposure of MMNK-1 cholangiocytes co-cultured with THP-1 macrophages, but not monocultured MMNK-1 cholangiocytes, to 1,2-DCP induced activation-induced cytidine deaminase (AID) expression, DNA damage and ROS production. The aim of this study was to identify relevant biological processes or target genes expressed in response to 1,2-DCP, using an in vitro system where cholangiocytes are co-cultured with macrophages. The co-cultured cells were exposed to 1,2-DCP at 0, 0.1 or 0.4 mM for 24 h, and then the cell lysates were assessed by transcriptome analysis. 1,2-DCP upregulated the expression of base excision repair genes in MMNK-1 cholangiocytes in the co-cultures, whereas it upregulated the expression of cell cycle-related genes in THP-1 macrophages. Activation of the base excision repair pathway might result from the previously observed DNA damage in MMNK-1 cholangiocytes co-cultured with THP-1 macrophages, although involvement of other mechanisms such as DNA replication, cell death or other types of DNA repair was not disproved. Cross talk interactions between cholangiocytes and macrophages leading to DNA damage in the cholangiocytes should be explored.

workers developed cholangiocarcinoma is estimated to range from 100 to 670 ppm 2 . The estimated range of 1,2-DCP exposure levels during the process of ink removal was reported to be 150-620 ppm 19 , which are comparable to occupational exposure levels to other organic solvents in poorly ventilated workplaces, which ranged from several hundreds to 1000 ppm 20,21 . To determine the equivalent 1,2-DCP concentrations to be used in our cell culture studies that match the above blood levels, we used the following assumptions; human blood: air partition coefficient of 10.7 1,22 11 . Based on these assumptions, we used 0.1 and 0.4 mM for 1,2-DCP concentrations in the present study, representing comparable levels to those found in workers exposed to 1,2-DCP.
PCA plots and volcano plots. Data normalization and differentially expressed genes (DEGs) identification between samples was performed by TCC package 24 . Then PCA analysis was performed with function "prcomp" in the "stats" package of R software 26 . Volcano plots of DEGs between groups were generated with "EnhancedVolcano" package 27 . Detection of co-expressed gene modules. Co-expressed gene modules among differentially expressed genes in 1,2-DCP-exposed MMNK-1 cholangiocytes and THP-1 macrophages were detected using Weighted Gene Co-expression Network Analysis (WGCNA) package 28 in R 3.5.1. Variance-stabilizing transformation of TCC-normalized count data was performed using the DESeq2 package 29 in R 3.5.1, and the transformed data were used as input in the WGCNA package. The power value used was 9 for MMNK-1 cholangiocytes and 10 for THP-1 macrophages, while the merge_thre value was 0.2, the threshold value for the output of co-expression interactions was 0.25, and other calculation settings were set to defaults in the WGCNA for both MMNK-1 cholangiocytes and THP-1 macrophages.
The genes in the gene module groups detected by WGCNA were further clustered into positively and negatively correlated gene groups by using pheatmap package in R 3.5.1.

Enrichment analyses.
Gene Ontology (GO) and the Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis of differentially expressed genes, were performed using the free online platform; the WEBbased gene set analysis Toolkit (WebGestalt), accessed at http:// www. webge stalt. org 30 . These were used to determine the over-representation (enrichment) analyses of the study with set parameters of a minimum of five genes and maximum of 2000 genes for a category, and False Discovery Rate (FDR) cut-off value of < 0.05, using Benjamini-Hochberg method for multiple test adjustment.
Selection of genes whose expression was dose-dependently changed by exposure to 1,2-DCP. Independently from cluster analysis and subsequent enrichment analysis, we used the Pearson correlation coefficient to assess the significance of dose-dependent changes in the expression of each gene following exposure to 1,2-DCP. In this analysis, the p-value of the difference in the expression level was adjusted using the Benjamini-Hochberg method 31  www.nature.com/scientificreports/ Expression of genes selected by hypothesis. In addition to the above comprehensive analysis, hypothesis-driven statistical tests were conducted. Since macrophages play diverse functions in the immune response to foreign substances and toxicants 32 and our previous studies showed 1,2-DCP-induced upregulation of AID, LDH cytotoxicity, DNA damage and ROS production in MMNK-1 cells only when they were co-cultured with THP-1-derived macrophages 11,13,14 suggesting involvement of intercellular signals, we conducted one-way analysis of variance (ANOVA) and post hoc Dunnett's multiple comparison for the expression levels of cytokines of TNF superfamily or interleukins, chemokines (CCL, CXCL, CL and CX3CL), cytokine/chemokine-related proteins and cytokine/chemokine receptors in MMNK-1 cholangiocytes or THP-1 macrophages after 24-h exposure to 1,2-DCP. We also conducted ANOVA and post hoc Dunnett's multiple comparison for expression levels of genes categorized for KEGG's base excision repair (BER), homologous recombination (HR) and non-homologous end joining (NHEJ) pathway. ANOVA and post hoc Dunnett's multiple comparison test were conducted using JMP Pro version 16.1.0 (SAS Institute Inc. Cary, NC).
Further analysis was conducted to evaluate the relevant processes and pathways involved in the transcriptomic profiles of co-cultures of MMNK-1 cholangiocytes and THP-1 macrophages exposed to 1,2-DCP for 24 h. For this purpose, overrepresentation analysis of gene ontology (GO) (biological process, cellular component and molecular function) and KEGG pathway terms, were employed. In MMNK-1 cells, significantly enriched (FDR < 0.05) term was detected only in the ME9 module ( Fig. 1f and Table 1). "Base excision repair" of KEGG pathway term was highly enriched in the ME9 module, showing overrepresentation of LIG1, PARP4, and POLD1 genes in the 1,2-DCP exposed group compared to the control group (Table 1).
In THP-1 macrophages, significantly enriched (FDR < 0.05) terms were detected only in the ME2 module (Fig. 2f, Tables 2 and 3). The cell cycle-related GO process/KEGG pathway terms were highly enriched in the ME2 module (Tables 2 and 3). The genes annotated to the KEGG pathway included BUBIB, CCNB2, CDC20, CDC45, CDK1, CDC7, MCM3, PLK1, and PTTG1, which were all upregulated in the 1,2-DCP group, compared to the control group (Table 3). The genes annotated to various GO terms for biological processes of the ME2 module in THP-1 macrophages included those related to the mitosis (Supplementary Table S1 Table S6) and NAD biosynthesis process (Supplementary Table S7). GO terms for cellular component  Table S8).

Hypothesis-driven analysis shows that inflammatory responses are induced by exposure to 1,2-DCP in THP-1 macrophages but not in MMNK-1 cholangiocytes.
Genes of the cytokines/ chemokines, cytokine-related proteins, cytokine/chemokine receptors, whose expression was significantly different between levels of exposure to 1,2-DCP, in co-cultured THP-1 macrophages included cytokines (TNFSF4, TNFAIP8L1, TNFAIP8L2-SCNM1 and receptor TNFRSF10A) and chemokines (CXCL2, CCL2, CCL7 and receptors CX3CR1, CCR6, CCRL2), but ANOVA showed significant difference between levels of exposure only in TNFAIP8L1, CCL2, CXCL2, CX3CR1 and CCR6 ( Table 6). None of interleukins or their receptors in THP-1 macrophages were significantly changed in expression. With regards to the cytokines, TNFAIP8L1 was downregulated in the 1,2-DCP group compared to the control group. Among the chemokines, CCL2 and receptors CX3CR1 were upregulated in the 1,2-DCP exposed group, compared to the control group, whereas CXCL2  Table S9). Regarding genes categorized for homologous recombination (HR) pathway, NBN and RPA1 were upregulated by exposure to 1,2-DCP, but no genes categorized for non-homologous end joining (NHEJ) pathway showed significant change in expression level between levels of exposure to 1,2-DCP.

Discussion
In this study, we investigated the transcriptomic profiles of co-cultures of MMNK-1 cholangiocytes and THP-1 macrophages after 24-h exposure to 1,2-DCP. We used the co-culture model of cholangiocytes and macrophages to mimic an inflammatory environment and exposed these groups of cells to 1,2-DCP to determine the transcriptional activities that occur under such an environment thereby identifying changes or processes occurring within cholangiocytes and macrophages leading to DNA damage, which is thought to play a pivotal role in carcinogenesis.
Immunohistochemical analysis of specimens of occupational cholangiocarcinoma showed high infiltration of inflammatory cells, even at sites of the bile duct in noncancerous hepatic tissues 3 . Importantly, Trush and Kensler reported increased toxicity of chemicals in the presence of inflammatory cells 33 . Furthermore, our group recently demonstrated the important role of macrophages in 1,2-DCP-induced cytotoxicity, reactive oxygen species production and DNA damage in cholangiocytes exposed to 1,2-DCP, which occurred only in the presence of macrophages 11,14 . As such, we sought to identify the transcriptional activities associated with the increased cytotoxic and genotoxic effects of 1,2-DCP, on co-cultured MMNK-1 cholangiocytes/THP-1 macrophages, to enhance our understanding of the molecular mechanisms involved in 1,2-DCP-induced cholangiocarcinoma.
In this study, we used 1,2-DCP concentration range comparable to the 1,2-DCP exposure levels experienced by the workers of the printing companies in Japan, who were diagnosed with occupational cholangiocarcinoma 2 , as described in the "Materials and methods" section.
KEGG pathway enrichment analysis showed base excision repair term was enriched in line with increase in 1,2-DCP level in the co-cultured MMNK-1 cholangiocytes (Table 1, Fig. 1e). Furthermore, it also showed  www.nature.com/scientificreports/ upregulation of all the genes annotated to base excision repair in the 1,2-DCP group compared to the control group (Table 1), consistent with our previous reports that DNA damage occurred in co-cultured MMNK-1 cholangiocytes following exposure to 1,2-DCP 11,13,14 .
Carcinogenesis occurs in three stages, namely: initiation, promotion, and progression. DNA damage has been established as the event that initiates carcinogenesis 34,35 . Faults in the DNA repair systems could also burden the cells with potential disadvantageous mutations 36 . More strand breaks could occur during the repair process, which could further enhance genomic instability or cell death 37 . It is therefore inferred that increased DNA damage can both enhance and compromise the survival of initiated cells when some damaged DNA escapes DNA repair or is left not fully repaired 34 .
In addition to DNA damage, initiation of cancers is enhanced in the presence of increased DNA damaging agent 34,38 . The upregulation of DNA repair genes suggests increase in DNA damage as 1,2-DCP concentration is increased, which could enhance mutation in the cells thereby increasing the resultant neoplasia. Moreover, the DNA damage in cholangiocytes co-cultured with THP-1 macrophages has been shown to be 1,2-DCP dose-dependent 14 . Immunohistological analysis of specimens obtained from the 1,2-DCP cholangiocarcinoma cases showed overexpression of γH2AX, a marker of DNA double-strand break, in the foci of BilIN, IPNB, invasive carcinoma, and non-neoplastic biliary epithelial cells, compared to specimen from control of common cholangiocarcinoma 7 .
The transcriptomic profiling of THP-1 macrophages co-cultured with MMNK-1 cholangiocytes exposed to 1,2-DCP indicated that enrichment of cell cycle related processes was associated with 1,2-DCP exposure (Tables 2  and 3). All the expression of the genes (BUBIB, CCNB2, CDC20, CDC45, CDK1, CDC7, MCM3, PLK1, and PTTG1) associated with the enriched terms were upregulated in the 1,2-DCP group, compared to the control group (Table 3). Because most of these genes are particularly engaged in ensuring the progression of the cell cycle from G1 to S and from G2 to M, and ensuring the proliferation of the cells 43,44 , exposure of co-cultures of THP-1 macrophages and MMNK-1 cholangiocytes to 1,2-DCP might induce the proliferation of the THP-1 macrophages. Because macrophages have the major role in the regulation of inflammatory responses and a subset of macrophages could locally proliferate 45 , accumulation of macrophages at the site of injury following exposure to 1,2-DCP possibly affects inflammatory responses, carcinogenesis, and tumor microenvironment. Table 4. Genes whose expression changed dose-dependently following exposure to 1,2-DCP in MMNK-1 cholangiocytes co-cultured with THP-1 macrophages. We tested the significance of Pearson correlation coefficient between the expression of each gene and 1,2-DCP level. The table lists the top or bottom five genes with the largest fold change at 0.4 mM in MMNK-1 cholangiocytes co-cultured with THP-1 macrophages when exposed to 1,2-DCP at 0, 0.1 or 0.4 mM for 24 h. Data of normalized values for gene expression are mean ± SD, n = 3. All p-values for gene expression were adjusted using Benjamini-Hochberg method and expressed as q-values. Fold change represents the value relative to the average of the control group (0 mM). www.nature.com/scientificreports/ Further analysis showed significant and dose-dependent changes in the expression of genes in the overall transcriptomic profiles of 1,2-DCP-exposed MMNK-1 cholangiocytes/THP-1 macrophages co-cultures (Tables 4  and 5). LIG1 (which was also found to be a component of the base excision repair pathway in KEGG analysis of MMNK-1 cholangiocytes) and FN1, which are implicated in diseases and cancer, were significantly correlated genes with the increase in 1,2-DCP level in the co-cultured MMNK-1 cells. DNA ligase 1 (LIG1) gene encodes a member of the ATP-dependent DNA ligase protein family, which plays a role in DNA replication, recombination, and repair pathways where it seals nicks in double stranded DNA 46 . Furthermore, previous studies demonstrated the engagement of LIG1 in various repair pathways, such as short-patch 47 or long-patch 48 base-excision repair, nucleotide excision repair 49 , mismatch repair 50 and non-homologous end-joining 51,52 . In pathological conditions, upregulation of LIG1 expression has been demonstrated in many human cancers 42 and mutations in LIG1 gene are associated with retarded joining of Okazaki fragments during DNA replication, hypersensitivity to a variety of DNA-damaging agents and aberrant DNA repair in human fibroblast strain (46BR) cells [53][54][55] . Fibronectin1(FN1) encodes a dimeric glycoprotein known to function in cell adhesion, cytoskeletal organization, migration, proliferation, and differentiation 56,57 . High FN1 levels have been associated with increased invasion and metastatic capability in lung and hepatic cancers 56,58 . It has also been reported to be a causative factor in the development of various pathological conditions, such as liver cirrhosis 59 . FN1 is also reported to stimulate the expression of various inflammatory factors in the tumor microenvironment, thereby highlighting the regulatory influence of this glycoprotein in major inflammatory cells 60,61 .
The results of hypothesis-driven gene expression analysis suggest the expression changes of TNF-α-induced proteins TNFAIP8L1, as well as chemokines CCL2 and receptors CX3CR1, CCR6 occurred in THP-1 macrophages ( Table 6) by exposure to 1,2-DCP, but not in MMNK-1 cholangiocytes. Our previous studies showed monocultured cholangiocytes exposed to 1,2-DCP showed no significant change in expression of γ-H2AX, suggesting the involvement of macrophages in the induction of increased DNA damage 14 . While the TNF-α related proteins remain strong candidates for extracellular signaling involved in DNA damage in cholangiocytes, further studies are needed to clarify their exact roles and cross talk between cholangiocytes and macrophages in 1,2-DCPinduced DNA damage in cholangiocytes.
Generally, activation of macrophages leads to the release of cytokines and chemokines, which contributes to crosstalk between the macrophages and their environment 62 . Since the primary function of cytokines is the regulation of immune and inflammatory responses of the host to the invading foreign substances or tissue injury, they play a vital role in ensuring the overall health of the host 63 . Repeated exposure to toxicants or xenobiotics induces Table 5. Genes whose expression changed dose-dependently following exposure of 1,2-DCP to THP-1 macrophages co-cultured with MMNK-1 cholangiocytes. We tested the significance of Pearson correlation coefficient between the expression of each gene and 1,2-DCP level. The table lists the top or bottom five genes with the largest fold change at 0.4 mM in THP-1 macrophages co-cultured with MMNK-1 cholangiocytes when exposed to 1,2-DCP at 0, 0.1 or 0.4 mM for 24 h. Data of normalized values for gene expression are mean ± SD, n = 3. All p-values for gene expression were adjusted using Benjamini-Hochberg method and expressed as q-values. Fold change represents the value relative to the average of the control group (0 mM). www.nature.com/scientificreports/ persistent production of cytokines and chemokines by macrophages, resulting in enhancement of inflammation, trying to control cellular stress and minimize cellular damage 64 . However, this could induce dysregulated cytokine and chemokine production, which could subsequently result in various pathological states and cancer 63,65 . We have demonstrated recently that 24-h exposure of THP-1 macrophages to 1,2-DCP results in upregulation of TNF-α, IL-1β and IL-6 protein expression 11,14 . In this study, the mRNA expression levels of TNF-α, IL-1β, and IL-6 were not significantly changed in the cell co-cultures. Previous studies have shown differences between expression changes in monocultured cells and when in co-culture with other cells 66 and differences between mRNA expression and protein expression 67 . However, the expression levels of certain tumor necrosis factor superfamilies of ligands (TNFSF) and receptors (TNFRSF) such as TNFSF4, and TNFRSF10A, and tumor necrosis factor-α-induced protein 8-like family (TIPE) such as TNFAIP8L1, TNFAIP8L2-SCNM1 were changed. TNFSF and TNFRSF are known to be expressed by or target immune cells such as macrophages and some non-immune cells through co-stimulatory and inhibitory pathways to induce the expression of wide range of actions including cellular differentiation, survival, and production of inflammatory cytokines and chemokines [68][69][70] . Active TNFSF ligand-receptor signaling pathways are associated with inflammatory disease and cancer 68,71 . More specifically, low expression of TNFSF4 mRNA was associated with worse prognosis in melanoma patients 71 . TIPE family has been described as regulators of immunity and tumorigenesis 72 . An expression analysis in humans showed that the TIPE family is dysregulated in cancer and inflammation and plays critical role in tumorigenesis and inflammatory responses 73 . TNFAIP8L1 has been demonstrated as an inducer of cell death 74 . Previous study showed the down regulation of TNFAIP8L1 in hepatocellular carcinoma (HCC) tissues which correlated with tumor pathogenic grade and patient survival 74 . In this study, interestingly 1,2-DCP exposed co-cultured macrophages showed a downregulation of TNFAIP8L1 in the 1,2-DCP exposed cells compared to the control group. Among genes categorized for base excision repair, only LIG1, PARP4 and POLD1 are differentially expressed in ME 9 but it should be noted that OGG1, which is categorized for base excision repair pathway, in ME7 is also upregulated by exposure to 1,2-DCP (Supplementary Table S9). Upregulation of base excision repair genes in cholangiocytes co-cultured with macrophages doesn't indicate that 1,2-DCP-induced DNA damage is exclusively limited to DNA base damage. Other than base excision repair pathway, LIG1 and POLD1 are also related to replication and other DNA repair pathways and PARP4 is also related to apoptosis or transportation as a vault protein. Our previous studies showed that exposure to 1,2-DCP increased DNA damage as assessed by alkaline comet assay, which detects both single strand breaks (SSBs) and double strand breaks (DSBs), and γH2AX expression, which detects DSBs, in cholangiocytes co-cultured with macrophages 11,14 . Table 6. ANOVA for expression levels of cytokines/chemokines-related genes or their receptors, which are selected by hypothesis, in THP-1 macrophages co-cultured with MMNK-1 cholangiocytes. Normalized values of expression level were compared between three groups of different 1,2-DCP concentration by one-way analysis of variance (ANOVA), being followed by post hoc Dunnett's multiple comparison with control (0 mM 1,2-DCP group). Data represents expression levels of cytokines/chemokines-related genes or their receptors from transcriptomic profiles of THP-1 macrophages co-cultured with MMNK-1 cholangiocytes when exposed to 1,2-DCP at 0, 0.1 or 0.4 mM for 24 h. Data of normalized values for gene expression are mean ± SD, n = 3. Fold change represents the value relative to the average of the control group (0 mM).

Genes involved in intercellular signal
Regulation Module eigengene p-value for ANOVA www.nature.com/scientificreports/ KEGG pathway enrichment analysis did not detect significant involvement of homologous recombination repair (HRR) or non-homologous end joining (NHEJ) pathway. When looking at each gene listed in any module eigengene, NBN and RPA1 in ME8 and POLD1 in ME9 categorized for HRR were upregulated by exposure to 1,2-DCP, but no genes categorized for NHEJ pathway in any module eigengene showed significant 1,2-DCP exposure-related change in expression level (Supplementary Table S9). Collectively, the study does not exclude possible involvement of HRR pathway with repair of DNA lesions.
On the other hand, the result did not detect 1,2-DCP-induced upregulation of AICDA expression. This might be due to the difference in the length of exposure between the present study and the previous study, as upregulation of AICDA was optimal after 9-h exposure to 1,2-DCP but fell down greatly after 12-h exposure to 1,2-DCP 13 , thus the present study does not disprove possible involvement of AICDA with 1,2-DCP-induced DNA damage in MMNK-1 cholangiocytes co-cultured with THP-1 macrophages. Interestingly a recent study shows that base excision repair is required for the processing of AID-induced lesions into DNA double strand breaks 75 . Further studies are needed to clarify the role of AID in 1,2-DCP-induced DNA damage in cholangiocytes.
The mechanism of how 1,2-DCP induces DNA damage has not been revealed. Our previous studies showed increase in reactive oxygen species (ROS) 14 , tail DNA% and tail moment in comet assay, or AID expression 13 in MMNK-1 cholangiocytes by co-culture with THP-1 macrophages, suggesting involvement of ROS or AID in DNA damage in cholangiocytes. Exposure to 1,2-DCP increased ROS level dose-dependently in MMNK-1 cholangiocytes co-cultured with THP-1 macrophages, but not in monocultured MMNK-1 cholangiocytes or THP-1 macrophages, suggesting ROS is produced by intrinsic mechanism in MMNK-1 cholangiocytes although it is activated by exposure to 1,2-DCP in the presence of macrophages 14 . Elevated ROS levels cause damage to DNA including abasic sites, single strand DNA breaks (SSBs), sugar moiety modifications, deaminated and adducted bases [76][77][78] . Oxidative base lesions such as highly mutagenic guanine derivative 7,8-dihydro-8-oxoguanine (8-oxoG) and the corresponding ring fragmented purine formamidopyrimidine derivative (FapyG) or abasic sites are predominantly repaired by base excision repair (BER) and to a lesser extent nucleotide excision repair 47,79,80 . Oxidative DNA lesions can lead to DNA double-strand break (DSB) formation which is originated from single strand break (SSB) during repair, excision of base, topoisomerase cleavage, DNA replication or transcription [81][82][83][84][85] . Upregulation of BER genes in the present study may be a response to ROS-induced DNA damage, although exact mechanism on how DSBs are generated is not clarified. Given the fact that exposure to 1,2-DCP increases the number of cholangiocytes with γH2AX-positive foci or the number of γH2AX-positive foci per nucleus of cholangiocytes, suggesting occurrence of DSB in cholangiocytes, in the presence of macrophages 11,14 , it is likely that various pathways of DNA damage repair may be involved. Studies using cells with pathway-specific geneknockouts are needed to fully understand how efficient but occasionally erroneous DNA damage repair occurs in cholangiocytes exposed to 1,2-DCP in the presence of macrophages.

Conclusions
The transcriptomic profiles of MMNK-1 cholangiocytes showed that the upregulation of base excision repair genes, and that such upregulation was 1,2-DCP-concentration dependent, indicating increased DNA damage in the cholangiocytes. The transcriptomic profiles of THP-1 macrophages, however; showed upregulation of cell cycle-related genes, indicating enhanced proliferation of macrophages. Upregulation of the base excision repair genes might be involved in the previously observed DNA damage in MMNK-1 cholangiocytes co-cultured with THP-1 macrophages, although involvement of other mechanisms such as DNA replication, cell death or other types of DNA repair was not disproved. Cross talk interactions between cholangiocytes and macrophages explaining the observed increase in DNA damage in the cholangiocytes should be explored further.

Data availability
Raw data, processed data and metadata of transcriptome analysis have been deposited in NCBI Gene Expression Omnibus (GEO; http:// www. ncbi. nlm. nih. gov/ geo; GSE 198858).